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Abstract 

Given a fluid equation with reduced Lagrangian I which is a functional of velocity u and 
advected density D given in Eulerian coordinates, we give a general method for semidiscretising 
the equations to give a canonical Hamiltonian system; this system may then be integrated in 
time using a symplectic integrator. The method is Lagrangian, with the variables being a set of 
Lagrangian particle positions and their associated momenta. The canonical equations obtained 
yield a discrete form of Euler-Poincare equations for I when projected onto the grid, with a new 
form of discrete calculus to represent the gradient and divergence operators. Practical symplectic 
time integrators are suggested for a large family of equations which include the shallow-water 
equations, the EP-Diff equations and the 3D compressible Euler equations, and we illustrate the 
technique by showing results from a numerical experiment for the EP-Diff equations. 
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1 Introduction 

There has been much recent interest in obtaining Hamiltonian methods for the various equations of 
motion for fluids (such as the shallow-water equations, and the 2D Euler equations). The Hamiltonian 
method programme consists of two stages: 

1. Take a Hamiltonian PDE with Hamiltonian H and structure operator J' so that the equation 
takes the form 

Discretise the Hamiltonian and the structure operator to obtain a finite dimensional Hamiltonian 
H and symplectic structure matrix J to give a finite dimensional Hamiltonian system 



2. Discretise the finite dimensional system in time using a symplectic integrator (for a general review 
of symplectic methods see (|HLW02|I or (fLROSj) ). 



The advantage of this approach is that the symplectic integrator will guarantee excellent long-time 
conservation properties IIBG94I [Hai94l IH L97l IRei99|l with the spatially-discrete Hamiltonian H satisfying 

\H{t)~H{S))\<ci5tP, \t\<C2e-'^'^\ 

where p is the order of the symplectic method in time and where ci, C2 and C3 are positive constants. 
Note that Hamiltonian methods are not the same as multi-symplectic integrators (MPS98 MPS VToTl 
[ReiOO M R03b MR03a ) where the Lagrangian is symmetrically discretised in time and space, and 
discrete variations are taken to obtain the numerical method, which then offers additional symmetry 
properties (OWW04). 

In general it has been very difficult to obtain discrete structure operators J which satisfy the Jacobi 
identity and progress has only been made for Hamiltonian PDEs with canonical structure (with a few 
exceptions e.g. fZei91')). In the context of fluid equations this means that the equations must be solved 
in Lagrangian coordinates where the symplectic structure is canonical. For the case of 2D incompressible 
flow this leads to methods where Lagrangian particles carry vorticity instead of mass, either as Dirac 
(5-functions (point-vortex methods (CKOO for a review)) or in "blobs" surrounding the particle (vortex- 
blob methods (Cho73 OSOl CR04bj). For the case of compressible flow the particles carry masses 
which are interpolated through basis functions to the entire domain to give the density field; this type 
of discretisation is called Smoothed-Particle Hydrodynamics (SPH) (GM77 Luc77 SalSS). The related 
Hamiltonian Particle-Mesh (HPM) method (FGR02 CFR03 FR 04_a CR0 4a) discretises a regularised 
form of the equations (with the regularisation operator discretised on an Eulerian grid) so that good 
long-time behaviour is obtained when the finite dimensional system is integrated with a symplectic 
integrator (this is described further in section 

The difficulty with Lagrangian coordinates is that the equations can become very complicated 
(especially when differential operators in Eulerian coordinates are involved). Analytically, the most 
tractable approach is to transform the Lagrangian L for the flow maps in Lagrangian coordinates to 
a reduced Lagrangian I for the velocity field u{x) in Eulerian coordinates, and then to use the Euler- 
Poincare (EP) equations ( H M R98a' ' H M R98E|l for u which are equivalent to the canonical Hamilton's 
equations for the Lagrangian variables. Writing the reduced Lagrangian which integrates over Eulerian 
coordinates means that the equations are much easier to handle. This greatly aids model derivation 
(pHoi02), analysis C MRSQOI IShkQ2|l . model reduction rHZ98l IHTiOSjl and averaging r Hol99|l : it 
is clear that it would be attractive to use the reduced Lagrangian to derive numerical methods too. 
However, in the case of general fluid equations this would mean losing the canonical structure operator 
because the Legendre transform in Eulerian coordinates results in a non-canonical Lie-Poisson equation. 

In this paper we give a general approach, using Lagrangian variables so that there is a canonical 
structure operator which is easy to discretise, whilst evaluating the Hamiltonian on an Eulerian grid. 
This approach fully develops a dual grid-particle approach (as an extension of the mixed form that HPM 
uses with the potential energy in Eulerian variables whilst the kinetic energy is in Lagrangian variables) 
The particle-mesh formulation allows quantities (such as momentum and density) which are given at 
the particle locations, to be interpolated to the Eulerian grid. The Hamiltonian is then approximated as 
a sum over the Eulerian grid points, and Theorem 13. 3l in this paper shows that the canonical equations 
obtained from the discrete Hamiltonian produce a set of discrete EP equations on the grid (although 
computationally the equations are treated in Lagrangian form). It then remains to apply a numerical 
time-stepping algorithm which preserves the canonical symplectic structure. 

The structure of the rest of the paper is as follows: section |2 provides some background on SPH 
and HPM to prepare the way for the general approach to obtaining semidiscretisations. This approach 
is set out, together with the statement and proof of Theorem 13.31 in sectional and suitable symplectic 
time-integrators are suggested in section |4l The next two sections show how to apply the approach in 
two examples: the shallow-water-a equations in section l5!T1 and the 2D EP-Diff equation in section 
15.21 Both of these examples are in 2-dimensional geometry although the approach is completely general 
and may be used to discretise n-dimensional PDEs. Section Ogives some numerical results for the 2D 
EP-Diff equation and section contains some brief final discussion. 
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2 Background 



HPM was originally introduced in (|FGR02|I as a semidiscretisation of the regularised layer-depth shallow- 
water (RLDSW) equations: 

Ut + u-Wu = -.gV(l - a'^Ay'^D, Dt + V ■ (Du) = 0. 

The success of this method is based on its "dual picture" nature, as the Lagrangian particle dynamics 
provides a canonical symplectic structure (which is easy to treat numerically using symplectic integra- 
tors) whilst the inverse modified Helmholtz operator takes its simplest form on the Eulerian grid. 

HPM can be viewed as a modification of SPH (GM771 ILuc77)l which is a Lagrangian particle 
method for compressible flow. In the SPH method for shallow-water, the layer-depth D is represented 
by a linear combination of radial basis functions which are centred on a finite set {^/3(i)}/3=i °f 
Lagrangian particles: 

1 " 

D{x,t) ^ —Y^DpcPix - Xpit)), 

13=1 

where (/> is a radial basis function with integral A5, and {i'/sj^^i are constant weights. Other fields 
f{x,t) (such as temperature, velocity etc.) may be interpolated from the values {fp{t)}''p^i at the 
particle locations to the entire domain via the product 

1 " 

13=1 

This interpolation may then be used to establish the continuity equation for h: 

dtD{x,t) = J^dtY,D0cl>{x-Xp{t)), 

p=i 



p=l 
1 " 

-AS ^ ^/'^Mi)Vx^(a; - Xpit)l 

13=1 
1 " 



0=1 

= ~V^-{Dv{x,t)), 

where Vx^ is the gradient with respect to Xp, is the gradient with respect to x, and v{x,t) is 
the interpolation of the particle velocities to the entire domain. 

In SPH methods, the forces on the particles are simply evaluated by interpolating the forces to 
the entire domain and then substituting the particle positions in the formula. These methods are thus 
mesh-free which makes them popular with astrophysicists who are able to model self-gravitating fluids 
(such as colliding stars and accretion discs) contained in an infinite vacuum simply by giving the basis 
functions finite support so that there is no density in the far-field. In particle-mesh methods, however, 
the aim is to meet the need for easy evaluation of Eulerian differential operators such as the inverse 
modified Helmholtz operator in the RLDSW equations. For a general mesh with m grid points {xkY^'^i, 
the layer-depth Dk at grid point k is given by 

1 " 

D{xk,t) ^—Y^D^M^p), 

(3=1 

where 

xPk{X) = c^{xk;X), 

and (f){x\Xi3) is some general basis function (such as a radial basis function or a Cartesian product of 
ID basis functions) centred on Xp that satisfies 

V^dp{x- X) - ~Vx(i}{x; X), J (p{x- 0) x = AS. 
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The key to the method is the restriction to basis functions (such as the Cartesian product of cubic 
B-splines ()dB78l for example)) which satisfy the partition of unity property: 

n 
13=1 

which means that any quantity / may be interpolated from the grid to the whole domain: 

m 

fix) ^^f{Xk)^pk{x). 
k=l 

This allows us to interpolate from the grid to the particles via 

m 

/^ = /(X^)=^/(a;fc)Vfc(X^). 

fc=i 

The gradient of the density can be obtained by taking the gradient of the basis functions evaluated 
at the particle positions: 

m 

V xD(x)\x=X 13 — 

fc=i 

m _j n 

fe=l 7=1 

The RLDSW equations were originally introduced after failed attempts to exploit the canonical 
Hamiltonian structure of SPH (applied to the SW equations) by performing the time integration using a 
symplectic method (with the aim of obtaining longtime approximate preservation of energy, momentum 
and adiabatic invariants). These attempts failed because computational efficiency necessitates the use 
of basis functions with compact support. When the flow is nearly incompressible (e.g. when the system 
is near to geostrophic balance in the rotating case) the basis functions supported local compressible 
oscillations which quickly destabilised the flow and led to equipartition of energy much too quickly. The 
smoothing operator applied to the layer depth means that the local basis functions become global and 
the small oscillations are filtered out in the momentum equation. The semi-discretised equations are: 

m m n 

Vp = g^Vx,^fc(X^)^(A-i)fe,— 5]Vfe(^7)^7, 

k=l 1=1 7=1 

= Vfj, (3=l,...,n, 

where A is the numerical approximation to the modified Helmholtz operator (obtained using e.g. Fourier 
series or finite elements). These equations have a canonical Hamiltonian structure with Hamiltonian 

n I 12 m n 

13=1 ^ k,l=l /3,7=1 

where the momentum is given by 

mp^DpVp, (3=l,...,n. 

These equations may be integrated using the symplectic Stormer-Verlet method because the Hamilto- 
nian splits into separate functions of {mj^^j^ and {Xj^^j^ (as in classical mechanics) (McL99). 

In this method, the potential energy is obtained by interpolating the height field to the grid and 
summing over the gridpoints, whilst the kinetic energy is evaluated in Lagrangian coordinates only. In 
the next section we shall extend the approach to the kinetic energy used here by evaluating the entire 
Hamiltonian on the Eulerian grid. 



4 



3 Canonical Hamiltonian particle-mesh semidiscretisations 



Throughout this section, and for the rest of this paper, we shall adopt the notation that a bar indicates 
a quantity stored at a particle location e.g. rn, D, etc., and a tilde indicates a quantity stored at a grid 
location e.g. m, D, u. Also, Greek letters will be used for indices denoting the particle numbering, 
e.g. momentum TtT/j at particle Xp, whilst Roman letters will be used for indices denoting the grid 
numbering, e.g. velocity Uk at grid point Xk. 

The programme for producing Hamiltonian particle-mesh semidiscretisations is as follows: 

1. Start with the reduced Lagrangian I for the EP equations given in Eulerian coordinates. This 
method applies to semi-direct product systems where the / is a functional of the velocity u and 
the density D (which satisfies the continuity equation). In this section we shall look at PDEs In 
two dimensions, but the extension to higher (and lower) dimensions is straightforward. 

2. To obtain the discrete Lagrangian L, take the continuous Lagrangian I which is written as an 
integral over P 

/ = J X{u,D)d^x, 
where A is a nonlinear operator, and replace it with a sum over gridpoints 

L = ^M«Afe {{ui}Zi,{Di}T=i), 

kl 

where M is the matrix representing the approximation to integration (often referred to in the 
finite-element literature as the mass matrix) and where is a function approximating A at the 
point Xk, with any differential operators being replaced by matrix operations on the grid in the 
manner consistent with the choice of M. 

3. Use the Legendre transform 

kl 

where 

^^^'^' = ^' 
to obtain the Hamiltonian written in Eulerian coordinates. 

4. Take n Lagrangian particles {-^/3}^=i together with m Eulerian grid points {xk}^^i on the 
domain V and a set of basis functions il>k{x) defined by 

ipk{x) = (l){xk;x), 

with the partition-of-unity property 

^iPkix) = 1 for all x€V. 

k 

and where (p satisfies 

Va,0(y;a;) = -Vj,0(y;a;), / (j){x;y) x = AS, 

Jv 

This property allows grid values f{xk) to be interpolated to the whole domain V: 

fi'^) = Yfixk)i^kix), 

k 

as well as gradients: 

Wf{x) = Y,fixk)^Mx)- 

k 

Here we avoid the issue of boundary conditions by stipulating that X> is a closed, compact manifold 
{e.g. a square with periodic boundaries). 
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5. Give the /3-th particle a momentum density and a constant layer-depth density Dp (which will 
be identified as a conserved mass in the particle system that is obtained after semidiscretisation). 
Interpolate the layer-depth density and momentum density to the grid using 

^ A/«m, = J2 ^MXp), E ^^fc'^' = E ^MXp). (2) 

10 I p 

The mass matrix is present because the product of the momentum m with the basis functions i/j 
should be interpreted as the approximation to the one-form density md^ x; inverting the mass 
matrix "strips off" the density x. 

6. The equations of motion for the particle positions {-X^/sj^^i and their associated local momenta 
{m^j/AS'l^^j^ are then the canonical Hamilton's equations for this discrete Hamiltonian JJl 
having applied the substitution in equations Q. 

Now we shall show that the equations of motion represent an approximation to the EP equations. 
First, we make a few definitions which set out our discrete calculus which is used to form the discrete 
EP equations. 

Definition 3.1. For a quantity f specified on the grid, define the grid-to-particle map by inter- 
polating to the whole of V, 

fix) = ^/feVfe(a;), 

ft; 

and evaluating at the particle positions 

k 

We shall use the following notation 

k 

Further define the grid-to-particle gradient map applied to f by interpolating to the whole of T), 

fix) = ^fktpkix), 

k 

taking a gradient, 

V/(a;)=^/fcV^fc(x), 
k 

and evaluating at the particle positions 

Vf{X0)=J2fk'^MXi3). 



We shall use the following notation 



V/ 



J^fk^MXp). 



Definition 3.2. For a density f specified on the particle locations, define the particle-to-grid map 
that interpolates to the grid by 

Y^Mu{f)i^Y.^MXp). 

I p 

Also for a function g specified on the grid, use the following notation for the interpolated product: 

J2Mki{7[g])i ^J^i^^ahMXp). 
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For a one-form density f specified on the particle locations, further define the particle-to-grid 
divergence map on the grid by 



AS' 



I /3 

This map is the adjoint to the grid-to-particle gradient map. 

These interpolation maps represent the whole key to the approach: they allow the dynamical equa- 
tions to be interpreted either on the grid or on the particle positions. In particular, the canonical 
symplectic structure is used to derive the particle dynamics (and so most of the focus computationally 
will be on the particle variables) whilst the Eulerian interpretation makes it possible to solve for the 
velocity given the momentum on the grid. The following theorem shows that the equations on the grid 
can be interpreted as a discrete EP equation: 

Theorem 3.3. On the grid the equations of motion obtained from this procedure take the following 
form 



-(m)fe + (V-(Nm)), + ([(V«r].m), = (d 



dD 



dL 



(3) 
(4) 
(5) 



which gives an approximation on the grid, using the discrete calculus defined above, to the EP 
equations: 



dtm + V • (um) + (Vm)^ • m = DV 



dtD + V ■ {uD) = 
m = 



0, 
5l_ 
5u 



Furthermore, at the particle locations the equations takes the form 



dD 



-1/3 



di' 



Xa = 



(6) 
(7) 



which gives an approximation at the particle locations to the EP equation in vector form: 

D_rn 

Wt'D 



Dm ,„ ,r m 
— - + (Vtxf •- 



D 
D 

Dt 



^5D' 



= {dt + u-V). 



Proof. The X equation obtained is 

Xf3 = ASVrn.H, 



X] I (^'"/s X] ^kirhi) ■ ilk + i'^mf^Uk) ■ X] ^kiirik j - Vwi^L. 
k \ I I / 



These three terms may be expanded out: 



X] ( X! ^kirhk \ Uk = X] ( X] iT^i'^kiXj) J • Uk 
k \ I ) k \ 1 ) 
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k 

k I kirn 

kl 



klm 

where we have made use of 

I 

This means that the second and third terms cancel and the X equation becomes 

Xfj = '^UkiJJk{Xp) = [u]p. 

k 

The m equation is 

= ^ {{'^ x„Mkirhi) ■ Uk + (Vx^Mfc) • Mkirfii) - Vx^i- 

kl 

Once again, expand these three terms 

kl k 

kl klm ^ ^ 



kl 



klm 



-Si:i;(«-)«vx.*(^rf. 



The m equation may then be written in the form 

^rfip = -[{Vuf]p ■ mp + Dp 

and dividing through by Dp gives equations H6I7|) . 

To obtain equation (PJ, take the time derivative of {rrT)k- 



dD 



-1/3 



dD 



(V-(m[«])),j, 

and dividing by M gives the result. Similarly, equation |(4j follows by taking the time derivative of 

□ 

This means that is possible to take any reduced Lagrangian l{u,D) and, following the procedure 
described above, obtain a finite-dimensional canonical Hamiltonian system with approximate Lagrangian 
L « Hor a set of particles, which may be interpreted on the grid as a discrete form of the EP equations 
resulting from /. If this system is integrated in time using a symplectic method, then the spatially- 
discretised Hamiltonian will be approximately preserved for long time-intervals. Suitable symplectic 
methods are discussed in the next section. 



4 Symplectic time integration methods 

in order to make a useful and practical Hamiltonian numerical method, it is necessary to choose the 
symplectic time-integrator carefully. In order to preserve the symplectic form, many integrators have 
to be implicit; often the expense of solving the resulting algebraic equations means that the symplectic 
integrator cannot compete with non-conservative schemes. In this section we suggest a practical first- 
order symplectic scheme for solving the semidiscretised system obtained using the method of the previous 
section. In order to do this, we restrict to a smaller set of problems where the Lagrangian takes the 
form 



J JC{D,u)-V{D)d'^x, 



where the nonlinear operator IC is interpreted as the kinetic energy and V is the potential energy; we 
shall further require that IC takes the form 

/C = -DuBu, 
2 

where S is a linear operator, so that m = B^^u. This set represents a very large number of fluid 
systems so is not too much of a restriction (in the examples we shall also consider the EPDiff equation 
which does not have the density factor). We shall proceed by discretising C: 



Vi(D) 



L = ^ Mki I ]^DiUi ■ ^(B-i);„,it,, 

kl \ m 

For a given time step size At, write 

« Xp{n^t), 

and similarly for m etc. To further compact the notation write 

im - E h^i^'p). (1)1 = J27pi^{x"^), 

k 

and adopt similar notation for the gradient operators. We shall also write Uk{[fn], D) as the solution 
to 

{m)k^V^,L{Ui{m),D),D), 
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i.e. the grid momentum-to-velocity map, so that 

ki = ^{m)i-Ui{{m),D). 
The (first-order) symplectic-Euler scheme then takes the form 



At 



-[(VJ7((m 



■n+1 \ 



■n+l 



-Df-, 



-D, 



2 (;d)" 
■ i9y — ■ 

V^((D)») 



J7 ((m"+i)", (7^)") 



At 



-[i7((m"+i)", £)")]" 



/3' 



which is implicit in m but not X. The algebraic system for m can be efficiently solved using the 
following iterative scheme: 



m 



-AtD 



13 



-AtDf, 



"13 



— 

oD 



where rnl'-' is the jth element of the sequence converging to mJl (with the sequence initialised with 



m^' — rri'p ). This results in inverting a sparse matrix (with number of non-zero elements propor- 
tional to the number of particles) each iteration (and for the case V = Q the matrix becomes block 
diagonal as the different particles decouple). 

A second-order symplectic discretisation can be obtained by using the second-order Lobatto llla-b 
partitioned Runge-Kutta method as discussed in (|Rei99|l with reference to adaptive symplectic methods 
for molecular dynamics. 



5 Examples 

5.1 Example: the shallow- water-a equations 

In this section we derive a Hamiltonian particle-mesh semidiscretisation for the shallow-water-a (SW-a) 
equations. These equations are obtained by decomposing Lagrangian particle trajectories into mean 
and fluctuating parts, averaging along those trajectories, and then applying the "frozen turbulence" 
Taylor assumption, before finally assuming isotropic, homogeneous turbulence statistics flHol99() . 
The model is derived from the reduced Lagrangian 



I 



D 



■\Vu\'^ - gD)A'^x, 



where a is the mean lengthscale for the fluctuations. The variational derivatives are 



^ = {D + V ■DV)u, 
du 



SD 2 ^' ' 



and substituting these derivatives into the EP equation 



dtm + {Vuf ■ 

ou 



V- ( u— 

ou 



SI 
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gives the equations 

dtm + u ■ Vm + ( Vit)^ • m + mV ■ u = DV [ - gD 



A+V-(Dm) = 0, 

m = {D-a^V-DV)u. (8) 

In these equations the momentum m is advected by the velocity u obtained by inverting equation 
The equations have a more familiar appearance in the form 

{dt+u- V)(l - a^D-^V ■ D)-^u + {Vuf ■ (1 - a^D-^V ■ D)-^u 

^^{-gD + \{\u\' + a^\Wu\^)y 
Dt + V- (uD) = 0, 

which become the shallow-water equations 

{dt + u- V)u = -gVD, Dt + \/ ■ (uD) = 0, 
when is set to zero, making use of the identity 

vi|n|^ = (Vm) ■ u. 

Following the procedure set out in section 13 we discretise the Lagrangian using a piecewise-linear 
Galerkin finite element discretisation (see (|Red93|l for example). The Galerkin expansion for a function 
/ given on grid points is 

/(a;) = ^iVfc(a;)/,, 

k 

where Nk is a continuous basis function which satisfies 

Under this scheme, the Lagrangian becomes 

i = ^ Qwfc • BkiUk - ^gDkMkiDi 

kl ^ 

where the Helmholtz matrix B takes the form 

B,, = / V DkMix) {N,{x)N,ix) + a'VN,{x) ■ VN,{x)) x, 
Jn ^ 

and the mass matrix M takes the form 

Af= / N^{x)N.i{x)d^ X. (9) 



Both B and M are sparse, well-conditioned, symmetric matrices so they are very quick to invert using 
a preconditioned conjugate gradients algorithm. 

Using the Legendre transform, the momentum rhk on the grid is 



dL 

I 

and the discrete Hamiltonian takes the form 



rhk = = BkiUk, 



^ = H\ ■ ^^'^^^ B-^)iMr, + ^gDkMkiDi 

kl \ m 
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The canonical equations derived from the Hamiltonian H are then 



Dp 



-[(V«f],.^-5[V^], 
Dp 



„ 1 (m) 1 , 
. 2 (i?) ^ , 



and the momentum equation on the grid take the form 

|(m), + (V • (N m)), + ( [{Vuf] ■ rn), = 



1 m . 
^ ^ 2 {D) 



We can also make a semidiscretisation for the Euler-a equations by removing the potential energy 
term and setting a constraint 

D{xk)^l, k=l,...,m, 

via a set of Lagrange multipliers which give the pressure at each grid-point. It is possible to construct 
symplectic integrators for constrained Hamiltonian systems of this type ( |Jay96| [LS94I IRei96|l and these 
methods have already been applied to enforcing an incompressibility constraint in (|CFR03|I . 



5.2 Example: the 2D EP-Diff equations 

The EP-Diff equations l|HM05|l are the generalisation to rt-dimensional space of the Camassa-Holm 
equation HCH93|I . and represent the geodesic equations of motion for the norm. They represent the 
prototype for the family of regularised equations such as Euler-a (Hol99' 'MRSOO), although they also 
have applications in computational anatomy and image processing ( MTY Q2. ■Mum98j) . 
The 2D EP-Diff equation is obtained from the reduced Lagrangian: 



2 



a 



:|V«|- 



where a is a constant lengthscale parameter, and so the equation of motion is 

dtm + u ■ Vm + (Vm)"^ ■ m + mV ■ m = 0, m — [1 — a^l\)u. 

Following the programme set out in section O with a Galerkin finite element approximation, the 
discrete Lagrangian takes the form 



i = ^ \uk ■ {Akiui). 



kl 



where the matrix A is the discrete inverse modified Helmholtz operator 

Am^ I Nk{x)Ni{x)+a^VNk{x)-VNi{x)A^x. 
Jn 

The discrete Hamiltonian in Eulerian coordinates on the grid is then: 

H = ^ ^rhk ■ ^{A^^)kmM^iThi, 

kl m 

where Mki is the finite element mass matrix given by equation 
The canonical equations for this Hamiltonian are then 



Dp 



= [U{{rn))],^[A-\m)]0, 
= -[iVU{{m))f]p-mp, 



and the equation for m is the EP-Diff equations represented in our discrete calculus. 
The symplectic Euler method with these equations is 



-[(VC/((r 



-n+l\n\\T]n 



)Y 



ip-^p 



+1 
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and the iterative scheme to solve for rn^'^^ is 

m;5+i^^+i = - Ai[(VJ7((m"+iJ)"))^];^ • m;+i'-''+\ (10) 

Using this scheme only requires the inversion of one 3x3 matrix for each particle per iteration which 
leads to a very efficient method. 

6 Numerical results 

To demonstrate the particle-mesh approach, we give results using the method for the 2D EP-Diff 
equation described section l5^ (numerical results and discussion of emergent behaviour for this system 
are given in ( HS03)). The equations are solved in a 27r x 2tt periodic square, with a = 0.3133. 
The discrete Lagrangian was obtained using piecewise-linear quadrilateral finite elements on a regular 
128x 128 grid and the particles positions were initialised 16 particles in each grid cell. The equations were 
integrated using Matlab using the HPM C-mex files (see http : //www . cwi .nl/projects/gi/HPM/) 
to perform the grid-particle operations. The finite element matrices B and M were stored as sparse 
matrices and inverted using conjugate gradients with incomplete Cholesky preconditioning, reaching a 
residual of less that 1 x 10~^ in less than 10 iterations. The iterative scheme l)10j) converged with a 
total error of less that 1 x 10^^ in less than 3 iterations for each timestep (with dt — 0.0204) during 
the experiment. 

The initial data had the velocity zero everywhere except for two thin lines which are set to collide. 
The plots of x-velocity in fiefures [Tl4l show the lines reconnecting after collision and the formation of 
thinner "peak-shaped" lines, as reported in (HS03). A plot of the discrete Hamiltonian during the 
numerical experiment is shown in figure|5l This plot illustrates that, as the timestepping method used 
is symplectic, the discrete Hamiltonian is conserved within 0{At) of the initial value at time i = for 
very long time-intervals. This is the real benefit of discretising the variational structure of the equations. 
Higher-order time-integrators will produce smaller errors in the Hamiltonian as A< ^ 0. 




Figure 1: Plot showing magnitude of velocity at i = 0. The velocity is distributed along two lines 
which are set to collide. 
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Figure 2: Plot showing magnitude of velocity a,t t — 0.2445. The two lines have stretched out and 
expanded before overlapping. 

7 Discussion and outlook 

We introduced a general framework for deriving canonical Hamiltonian semi-discretised equations from 
the reduced Lagrangian in Eulerian coordinates, based on the particle-mesh approach. This framework 
has several advantages: 

• The use of Lagrangian particles in the discretisation means that the equations are always canonical. 
This means that the discrete structure operator is guaranteed to satisfy the Jacobi identity and so 
a symplectic time integrator can be used, leading to long-time preservation of the Hamiltonian. 
The canonical structure also makes for easier discussion of adiabatic invariants in the discrete 
system (|Cot04|l . 

• The Hamiltonian is given in terms of Eulerian coordinates. This means it is easy to add on 
extra physics, make approximations, etc. Many models, such as those discussed in section |^ 
involve differential operators given in Eulerian coordinates which become rather complicated and 
"entangled" in the Lagrangian variables. This approach avoids this problem. 

• This approach allows the momentum densities to become (5-functions in the limit, whilst the 
velocities on the grid can be represented using finite elements which means that weak solutions 
can be discussed, for example in the EP-Diff equations we seek solutions in (as illustrated in 
the numerical example). 

• The symplectic Euler method can be efficiently applied to integrate the spatially-discrete equations 
in space. It should also be possible to obtain higher-order symplectic methods for greater accuracy. 

• The numerical experiment demonstrates that, at least for the EP-Diff equations, this numerical 
method is computationally feasible and produces well-behaved numerical solutions. 

Future work may take several different directions, including: 

• Development of higher-order symplectic integrators for these models in order to produce practical 
methods. 

• Application to other fluid models: the two-layer Green-Naghdi equations f CC99|l . the quasi- 
geostrophic equations (.HZ98,1 . rotating shallow-water-a, Euler-Boussinesq-a in 3D (|Hol99|l etc. 
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Figure 3: Plot showing magnitude of velocity at t = 0.4074. The lines have reconnected after 
crossing. 

• Comparison of methods with existing schemes for simple models such as the EP-Diff equations 
in 1, 2 and 3 dimensions. 

• Investigation of the effects of smoothing in the shallow-water-a equations. Does the smoothing 
in the velocity field keep the solution well-behaved over long time intervals? 

• Use of normal-form theory, following the work in ^Cot04|l . to investigate geostrophic balance in 
the numerical scheme for the rotating shallow-water-a equations, and other rotating models. 

• Comparison of the solutions of shallow-water-a obtained using the method given in this paper 
with solutions of RLDSW obtained using the HPM method. 

• Inclusion of boundary conditions: the approach, as described here, only works on a closed manifold 
(such as the 2-torus i.e. a square with periodic boundary conditions). 

• Obtain an extension to semi-direct product systems with other types of advected quantities (other 
than densities) such as scalars, vectors, and higher-order tensors. 

Acknowledgements: Thanks to Sebastian Reich, Darryl Holm and Joel Fine for their useful discussions 
during the writing of this paper, to Jason Frank for the use of his HPM C-mex code, to Matthew 
Piggott for advice on programming the finite element method, to the Daniels family of Ohio for their 
generous hospitality during the early stages of this work, and to NERC for funding under grant number 
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